package project2;

/*********************************************************************
 * Get an optimal global alignment in linear space using Hirschberg's trick
*********************************************************************/

public class Hirschberg_global {


     /*********************************************************************
     * score_matrix[][] : matrix of 4x4
     * gap_cost : the cost for a gap in the alignment
     *********************************************************************/

     int[][] score_matrix;
     int gap_cost;

    Hirschberg_global(int[][] m, int g) {
        score_matrix = m;
        gap_cost = g;
    }


    /*********************************************************************
    * used to compute score_matrix[value(s1.charAt(i-1))][value(s2.charAt(j-1))]
    *********************************************************************/

    int value(char c){
        if(c=='a') return 0;
        if(c=='c') return 1;
        if(c=='g') return 2;
        if(c=='t') return 3;
        return -1;
    }


    /*********************************************************************
    * global(s1, s2) finds an optimal global alignment in quadratic space
    *********************************************************************/

    String[] global(String s1, String s2) {
        int x = s1.length();
        int y = s2.length();
        String[] result = new String[2];
        int[][] cost = new int[x+1][y+1];
        String op = ""; //the string I use for keeping track of the operations
        for(int i=0; i<=x; i++) {
            for(int j=0; j<=y; j++) {
                int v1, v2, v3, v4;
                v1 = v2 = v3 = v4 = 50000;
                if((i>0)&&(j>0))
                    v1 = cost[i-1][j-1] + score_matrix[value(s1.charAt(i-1))][value(s2.charAt(j-1))];
                if((i>0)&&(j>=0))
                    v2 = cost[i-1][j] + gap_cost;
                if((i>=0)&&(j>0))
                    v3 = cost[i][j-1] + gap_cost;
                if((i==0)&&(j==0))
                    v4 = 0;
                cost[i][j] = java.lang.Math.min(java.lang.Math.min(v1, v2),java.lang.Math.min(v3, v4));
            }
        }
        //compute the optimal operations
        int i,j;
        i = x; j = y;
        while((i!=0)||(j!=0)){
            if(i == 0) {
                j--;
                op = op.concat("I");
            }
            else
                if(j == 0) {
                    i--;
                    op = op.concat("D");
                }
                else {
                int v1, v2, v3;
                v1 = cost[i-1][j-1] + score_matrix[value(s1.charAt(i-1))][value(s2.charAt(j-1))];
                v2 = cost[i-1][j] + gap_cost;
                v3 = cost[i][j-1] + gap_cost;
                int min = java.lang.Math.min(java.lang.Math.min(v1,v2),v3);
                if(v1==min) {
                    op = op.concat("M");
                    i--;
                    j--;
                }
                else
                    if(v2==min) {
                        op = op.concat("D");
                        i--;
                    }
                    else
                        if(v3==min) {
                            op = op.concat("I");
                            j--;
                        }
                }
        }
        result[0] = Integer.toString(cost[x][y]);
        result[1] = op;
        return result;
    }


    /*********************************************************************
    * finds the middle edge for Hirschberg's trick
    *********************************************************************/

    int[] find_edge(String seq1, String seq2, int l1,int l2,int r1, int r2,int h) {
        int[] edge = new int[2];
        int[] cost,trace,auxc,auxt,type,auxy;
        cost = new int[r2-r1+1];
        trace = new int[r2-r1+1];
        auxc = new int[r2-r1+1];
        auxt = new int[r2-r1+1];
        type = new int[r2-r1+1];
        auxy = new int[r2-r1+1];
        for(int i=0; i<=r2-r1; i++)
            cost[i] = trace[i] = auxc[i] = auxt[i] = type[i] = auxy[i] = 0;
        for(int i=0; i<=l2-l1; i++) {
            for(int j=0; j<=r2-r1; j++) {
                //compute the value of the cost
                int v1, v2, v3, v4;
                v1 = v2 = v3 = v4 = 50000;
                if((i>0)&&(j>0))
                    v1 = auxc[j-1] + score_matrix[value(seq1.charAt(i+l1-1))][value(seq2.charAt(j+r1-1))];
                if((i>0)&&(j>=0))
                    v2 = auxc[j] + gap_cost;
                if((i>=0)&&(j>0))
                    v3 = cost[j-1] + gap_cost;
                if((i==0)&&(j==0))
                    v4 = 0;
                cost[j] = java.lang.Math.min(java.lang.Math.min(v1, v2),java.lang.Math.min(v3, v4));
                //compute the traceback
                if(i<(h-l1)) {
                    trace[j] = -1;
                    type[j] = -1;
                }
                else
                    if(i==(h-l1)) {
                        if(cost[j] == v1){
                            trace[j] = j+r1;
                            type[j] = 1; //diagonal
                        }
                        else
                            if(cost[j] == v2) {
                                trace[j] = j+r1;
                                type[j] = 0; //vertical
                            }
                            else
                                if(cost[j] == v3) {
                                    trace[j] = trace[j-1];
                                    type[j] = trace[j-1];
                                }
                    }
                    else
                        if(cost[j] == v1) {
                            trace[j] = auxt[j-1];
                            type[j] = auxy[j-1];
                        }
                        else
                            if(cost[j] == v2) {
                                trace[j] = auxt[j];
                                type[j] = auxy[j];
                            }
                            else
                                if(cost[j] == v3) {
                                    trace[j] = trace[j-1];
                                    type[j] = type[j-1];
                                }
            }
            for(int j=0; j<=r2-r1; j++) {
                auxc[j] = cost[j];
                auxt[j] = trace[j];
                auxy[j] = type[j];
            }
        }
        edge[0] = trace[r2-r1];
        edge[1] = type[r2-r1];
        return edge;
    }


    /*********************************************************************
    * apply Hirschberg's trick to find an optimal alignment of
    * seq1[l1...l2] and seq2[r1...r2]
    *********************************************************************/

    public String globalH(String seq1, String seq2, int l1, int l2, int r1, int r2) {
        String op = "";
        //first treat the base cases
        if(l2-l1<3||r2-r1<3) {
            if(l1>0)
                l1--;
            if(r1>0)
                r1--;
            String[] result = global(seq1.substring(l1,l2),seq2.substring(r1,r2));
            op = result[1];
        }
        //deal with the general case
        else {
            int h = l1+(l2-l1)/2;
            int[] edge;
            edge = find_edge(seq1,seq2,l1,l2,r1,r2,h);
            if(edge[1] == 0) { //vertical
                String op1 = globalH(seq1,seq2,l1,h-1,r1,edge[0]);
                String op2 = globalH(seq1,seq2,h+1,l2,edge[0]+1,r2);
                op = op2.concat("D");
                op = op.concat(op1);
            }
            else { //diagonal
                String op1 = globalH(seq1,seq2,l1,h-1,r1,edge[0]-1);
                String op2 = globalH(seq1,seq2,h+1,l2,edge[0]+1,r2);
                op = op2.concat("M");
                op = op.concat(op1);
            }
        }
        return op;
    }


    /*********************************************************************
     * apply Hirschberg's trick to find an optimal alignment of seq1 and seq2
    *********************************************************************/

    public String Hirschberg(String seq1, String seq2) {
        String aux = globalH(seq1,seq2,0,seq1.length(),0,seq2.length());
        String op = new String("");
        //I get the opertaions in the reversed order so I have to reverse the obtained string
        for(int i=aux.length()-1; i>=0; i--)
            op = op.concat(Character.toString(aux.charAt(i)));
        return op;
    }

}
